%% atm_test.m  4/3/2014  Parker MacCready
%
% development code to read in and process one day of WRF forecast
%
% I am trying a clearer workflow structure as an experiment that could work
% for any of these "reformatting" jobs that come up so often.
%
% 1. Create vectors of output times that map to vectors of input files
% 2. For each time:
%   - get the input file
%   - read in all the variables
%   - process as needed (rotate, rescale, make relative humidity)
%   - regrid to the output grid (here just horizontal)
% 3. Write each variable to its place in the netcdf output file

clear; addpath('../../../alpha'); Ldir = Lstart;

% define the time/filename vectors

% the forecast day
yr = 2014; mo = 3; dy = 27;
td_now = datenum(yr,mo,dy);
td_prev = td_now - 1;
td_next = td_now + 1;

% create the output time vector
% assume we want 37 hourly fields, hours 00 through 36
hr_vec = 0:36;
dt_out = datenum(yr,mo,dy,hr_vec,0,0);
ts_out = (dt_out - datenum(1970,1,1))*86400; % seconds since 1/1/1970

%% create the corresponding input files
%
% for now we will practice by taking the first 12 hours of the 00 forecast
% and then 24 hours of the 12 forecast
%
% and for now let's just use the d2 (12 km) grid

yrs = datestr(td_now,'yyyy');
mos = datestr(td_now,'mm');
dys = datestr(td_now,'dd');
% input directory(s)
indir00 = [Ldir.parent,'LiveOcean_notes/wrf/',yrs,mos,dys,'00/'];
indir12 = [Ldir.parent,'LiveOcean_notes/wrf/',yrs,mos,dys,'12/'];

% make the infile list
for tt = 1:length(hr_vec)
    hr = hr_vec(tt);
    if hr <= 12
        hrs = ['00',num2str(hr)]; hrs = hrs(end-1:end);
        infile_list_d2{tt} = [indir00,'wrfout.ocean_d2.',yrs,mos,dys, ...
            '00.f',hrs,'.0000'];
    else
        hrs = ['00',num2str(hr - 12)]; hrs = hrs(end-1:end);
        infile_list_d2{tt} = [indir12,'wrfout.ocean_d2.',yrs,mos,dys, ...
            '12.f',hrs,'.0000'];
    end
    forecast_hour{tt} = hrs; % used for rain calculation
end
% nice!

%% variable name lists 

invar_list = {'PSFC','RAINC','RAINNC','SWDOWN','GLW', ...
    'T2','Q2','U10','V10'};


%% place for output

fname = ['f',yrs,'.',mos,'.',dys];
outdir0 = [Ldir.out,Ldir.gtag,'/',fname,'/'];
if ~exist(outdir0,'dir'); mkdir(outdir0); end;
outdir = [outdir0,'atm/'];
if ~exist(outdir,'dir'); mkdir(outdir); end;

%% grids

% the ROMS grid
gdir = [Ldir.res,Ldir.gtag,'/'];
fng = [gdir,'grid.nc'];
lon = nc_varget(fng,'lon_rho');
lat = nc_varget(fng,'lat_rho');

% the WRF grids (2 for d2, and 3 for d3)
lon2 = double(nc_varget(infile_list_d2{1},'XLONG'));
lat2 = double(nc_varget(infile_list_d2{1},'XLAT'));

% make blank results arrays
[NR,NC] = size(lon);
[NR2,NC2] = size(lon2);
NT = length(hr_vec);
nmat = NaN * ones(NT,NR,NC); % sized for output
nmat2 = NaN * ones(NT,NR2,NC2); % sized for input

%% read input files into separate variable files

% start with all variables at one time, and end up with
% one variable at all times

for tt = 1:NT
    fn2 = infile_list_d2{tt};
    for vv = 1:length(invar_list) 
        VR = invar_list{vv};
        if tt == 1; eval([VR,'2 = nmat2;']); end;
        eval([VR,'2(tt,:,:) = nc_varget(fn2,VR);']);
    end
end

%% processing

% rain
%
% WRF reports the accumulation since the start of the
% forecast, so we can estimate the precipitation RATE by
% dividing by the time since the start of the forecast.
% The 0.1 gets us from mm/sec to cm/sec

RAIN2 = nmat2; % precipitation rate [cm s-1]
for tt = 1:NT
    dts =3600 * str2num(forecast_hour{tt});
    if dts == 0
        dts = 3600;
        RAIN2(tt,:,:) = 0.1 * (RAINC2(tt+1,:,:) + RAINNC2(tt+1,:,:))/dts;
    else
        RAIN2(tt,:,:) = 0.1 * (RAINC2(tt,:,:) + RAINNC2(tt,:,:))/dts;
    end
end

% calculate relative humidity
%
T22 = T22 - 273.15; % convert K to C
RH2 = Z_wmo_RH(PSFC2/100,T22,Q22); % relative humidity [%]

% rotate velocity to E-W and N-S
%
theta = NaN * lon2;
for rr = 1:NR2
    [dist,theta(rr,1:end-1)] = sw_dist(lat2(rr,:),lon2(rr,:),'km');
end
theta(:,end) = theta(:,end-1); % pad the end
alpha = - pi*theta/180; % alpha is the angle (radians) to rotate the
% coordinate system to go from the WRF grid to East-North
sa = sin(alpha); ca = cos(alpha);
ssa = nmat2; cca = nmat2;
for tt = 1:NT
    ssa(tt,:,:) = sa;
    cca(tt,:,:) = ca;
end
UR102 = U102.*cca + V102.*ssa;
VR102 = V102.*cca - U102.*ssa;

%% regrid to the model grid

invar_list = {'PSFC','RAIN','SWDOWN','GLW', ...
    'T2','RH','UR10','VR10'};

outvar_list = {'Pair','rain','swrad','lwrad_down', ...
    'Tair','Qair','Uwind','Vwind'};

for tt = 1:NT
    for vv = 1:length(invar_list) 
        VR = invar_list{vv};
        vr = outvar_list{vv};
        if tt == 1; eval([vr,' = nmat;']); end;
        eval(['OO = squeeze(',VR,'2(tt,:,:));']);
        F = TriScatteredInterp(lon2(:),lat2(:),OO(:));
        eval([vr,'(tt,:,:) = F(lon,lat);']);
    end
end

%% write out to NetCDF for ROMS



